A nonlinear theory of non-stationary low Mach number channel flows of freely cooling 

nearly elastic granular gases 
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We employ hydrodynamic equations to investigate non-stationary channel flows of freely cooling 
dilute gases of hard and smooth spheres with nearly elastic particle collisions. This work focuses on 
the regime where the sound travel time through the channel is much shorter than the characteristic 
cooling time of the gas. As a result, the gas pressure rapidly becomes almost homogeneous, while 
the typical Mach number of the flow drops well below unity. Eliminating the acoustic modes and 
employing Lagrangian coordinates, we reduce the hydrodynamic equations to a single nonlinear and 
nonlocal equation of a reaction-diffusion type. This equation describes a broad class of channel flows 
and, in particular, can follow the development of the clustering instability from a weakly perturbed 
homogeneous cooling state to strongly nonlinear states. If the heat diffusion is neglected, the 
reduced equation becomes exactly soluble, and the solution develops a finite-time density blowup. 
The blowup has the same local features at singularity as those exhibited by the recently found 
family of exact solutions of the full set of ideal hydrodynamic equations (Fouxon et al. 2007). The 
heat diffusion, however, always becomes important near the attempted singularity. It arrests the 
density blowup and brings about novel inhomogeneous cooling states (ICSs) of the gas, where the 
pressure continues to decay with time, while the density profile becomes time-independent. The ICSs 
represent exact solutions of the full set of granular hydrodynamic equations. Both the density profile 
of an ICS, and the characteristic relaxation time towards it are determined by a single dimensionless 
parameter C that describes the relative role of the inelastic energy loss and heat diffusion. At C 1 
the intermediate cooling dynamics proceeds as a competition between "holes" : low-density regions 
of the gas. This competition resembles Ostwald ripening (only one hole survives at the end), and 
we report a particular regime where the "hole ripening" statistics exhibits a simple dynamic scaling 
behavior. 

PACS numbers: 45.70.Qj, 47.20.Ky 



I. INTRODUCTION 

Clustering of matter is a spectacular example of struc- 
ture formation in nature. A fascinating example of clus- 
tering is provided by granular gases: gases of macroscopic 
particles that lose kinetic energy in collisions. Granular 
gas is a low-density limit of granular flows 1, 2] . The sim- 
plest version of the granular gas model assumes a dilute 
assembly of identical smooth hard spheres (with diame- 
ter a and unit mass) who lose energy at binary collisions 
in such a way that the normal component of the rela- 
tive velocity of the colliding particles gets reduced by 
a constant factor < r < 1 (the coefficient of normal 
restitution) upon each collision. Granular gases exhibit 
various pattern forming instabilities, including the shear- 



instability causes the generation of a macroscopic flow, 
both solcnoidal and potential, and formation of dense 
clusters of particles. 

A natural theoretical description of macroscopic gran- 
ular flows is provided by the Navier-Stokes granular hy- 
drodynamics l2| . Although the criteria of its validity 
are quite restrictive, see below, granular hydrodynamics 
has a great predictive power, sometimes going far beyond 
the formal limits of applicability 2]. Recently, granu- 
lar hydrodynamics has been applied to a variety of non- 
stationary flows of granular gases [H, HE EE (3- 



Non-stationary flows provide sharp tests to continuum 
models of granular flows, especially when the time- 
dependent solutions of the continuum equations tend to 
develop finite-time singularities. Examples are provided 
by the recently predicted finite-time blowup of the gas 
density in freely cooling granular gases: at zero gravity 
[H, 03, (3 (as described by ideal granular hydrodynamic 
equations), and at finite gravity (even in the framework 
of non-ideal granular hydrodynamic equations) [l6j . 

We will assume in this paper that particle collisions 
are almost elastic, the local gas density (that we denote 
by p) is much smaller than the close-packing density, and 
the Knudsen number is very small: 

l-r<l, pcr d <l, and J/ ree /L«l. (1) 

Here d > 1 is the dimension of space, I free is the mean 
free path of the particles, and L is the characteristic 
length scale of the hydrodynamic fields. Under these as- 
sumptions (the second and third ones need to be verified 
a posteriori) the Navier-Stokes hydrodynamics provides 
a quantitatively accurate leading-order theory 0, ■ It 
was shown 0, by using hydrodynamic equations that, 
for sufficiently large systems, the homogeneous cooling 
state (HCS) of the granular gas becomes unstable with 
respect to small perturbations. There are two linearly un- 
stable modes. The shear mode corresponds to the devel- 
opment of a macroscopic solenoidal flow, while the clus- 
tering mode corresponds to the development of a macro- 
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scopic potential flow that brings about formation of clus- 
ters of particles. 

A consistent nonlinear hydrodynamic theory of the 
clustering instability has not been available for quite 
a long time. Solving strongly nonlinear hydrodynamic 
equations is hard (even numerically), and one looks 
for additional simplifications. Following Refs. [H, EH, 
UtI [l8j . we will assume throughout this paper that the 
macroscopic flow (but not microscopic motion of the par- 
ticles!) is one-dimensional (Id). This assumption is nat- 
ural in the geometry of a narrow channel with perfectly 
elastic side walls that we adopt here. In a narrow chan- 
nel both the clustering mode in the transverse directions, 
and the shear mode are suppressed (see Refs. EH, EH for 
detail). As a result, the macroscopic flow can depend 
only on the coordinate along the channel and time, and 
we can focus on the development of the pure clustering 
mode as it enters a strongly nonlinear regime. 

Efrati et al. [12] investigated numerically the long- 
wavelength limit of such a quasi-ld clustering instabil- 
ity. In this limit the inelastic energy loss of the gas is 
the fastest process, so the gas pressure rapidly drops to a 
very small value. The further dynamics becomes (almost) 
purely inertial which (almost) brings about a finite-time 
blow-up of the velocity gradient and, therefore, of the 
density . The signatures of this finite-time singular- 
ity were indeed observed in the numerical solution of the 
hydrodynamic equations [l^] until the growing gas den- 
sity became so high that the numerical scheme lost ac- 
curacy. The numerical results of Ref. I12fl were tested 
in molecular dynamics (MD) simulations [13}. The MD 
simulations supported the free-flow blow-up scenario un- 
til the time when the gas density approached the hexag- 
onal close-packing value, and the further density growth 
stopped. 

Recently, Fouxon et al. [U EH analyzed, analytically 
and numerically, the one dimensional flow in the frame- 
work of equations of ideal hydrodynamics (that is, ne- 
glecting the heat diffusion and viscosity effects). They 
derived a family of exact solutions to these equations, 
with and without shocks, for which an initially smooth 
flow develops a finite-time density blowup. Close to the 
blow-up time t c , the maximum density exhibits a power 
law behavior ~ (t c — t)~ 2 . The velocity gradient blows 
up as ~ — (t c — i) -1 , whereas the velocity itself remains 
continuous and develops a cusp, rather than a shock dis- 
continuity, at the singularity. The gas temperature van- 
ishes at the singularity, but the pressure remains finite. 
Extensive numerical simulations with the ideal hydrody- 
namic equations showed that the singularity exhibited by 
the exact solutions is universal, as it develops for generic 
initial conditions. Very recently, the existence of the at- 
tempted blowup regime has been proved in molecular dy- 
namic simulations of a gas of nearly elastically colliding 
hard disks in a channel geometry [2(J] . The results of Refs. 
[I?], EH also imply that, for long wavelength initial con- 
ditions, the free flow re gim e may not hold all the way to 
the density blowup [l7lll8j . Very close to the attempted 



free-flow singularity, compressional heating starts to act. 
As a result, the gas pressure again becomes important 
and changes the local blowup properties. 

A crucial feature of the finite-time singularity of the 
ideal hydrodynamic equations is that it obeys an iso- 
baric scenario: the (finite) gas pressure becomes uniform 
in space in a close vicinity of the developing singularity 
fi"H . This hints at the possibility of an additional simplifi- 
cation of the problem. Indeed, an (almost) homogeneous 
pressure in a gas implies a low Mach number flow, when 
the inertial terms in the momentum equation are small 
compared to the pressure gradient term. This regime ap- 
pears when the sound travel time through the system is 
very short compared with other time scales of the prob- 
lem, and one is interested in the dynamics of the system 
at the long time scales [HI, [H, Ha, H3, IH, [H, [23, HI] . 
In particular, this regime appears naturally in the linear 
theory of the clustering instability of the HCS for inter- 
mediate wavelengths of the perturbations, see below. It is 
this (almost) spatially independent pressure regime that 
we will be considering in the present work. 

The remainder of the paper is organized as follows. In 
Section [TI] we start with a full set of equations of granular 
hydrodynamic for a dilute granular flow in a channel and 
reduce them, for sufficiently short channels, to the low 
Mach number flow equations. In Section lllll we employ 
Lagrangian coordinates which enable us to exactly reduce 
the low Mach number flow equations to a single nonlin- 
ear and nonlocal equation, of a reaction-diffusion type, 
for the square root of the inverse gas density. The new 
equation is tested in Section ITvl on two simple problems: 
the HCS and the linear theory of clustering instability 
in short channels. In Section [V] we show that, when the 
heat diffusion is neglected, the new equation becomes 
exactly soluble, and the solution develops a finite-time 
density blowup with the same universal features at singu- 
larity as those exhibited by the family of exact solutions 
of the full set of ideal granular hydrodynamic equations 
[l7l[l8j. Section IVTl presents an analytical and numerical 
analysis that shows that the heat diffusion term, no mat- 
ter how small in the beginning, becomes important near 
the attempted density blowup. As a result, the density 
blowup is arrested, and a novel, inhomogeneous cooling 
state (ICS) of the gas emerges, with a time-independent 
inhomogeneous density profile. Importantly, the ICSs 
represent exact solutions of hydrodynamic equations. A 
limiting form of the novel cooling state is what we call the 
"hole" , and we investigate its properties and the relax- 
ation dynamics towards it. For sufficiently long channels 
(other parameters being fixed) the cooling dynamics of 
the system takes the form of a competition between, and 
"ripening" of, holes. Therefore, in Section fVIII we inves- 
tigate the dynamics and statistics of this competition. In 
Section lVIIII we summarize our results and put them into 
a perspective. 
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II. GRANULAR HYDRODYNAMICS AND A 
LOW MACH NUMBER FLOW 

For flows depending on a single spatial coordinate x 
and time t the granular hydrodynamic equations can be 
written as follows: 



dp d(pv) 
dt dx 
' dv dv 

P ' 
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dv\ 
r ^~dx) ' 



«8_ f^dT 
p dx \ dx 
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1)VT fdvY 
[dx] 



(4) 



Here 7 is the adiabatic index of the gas (7 = 2 and 5/3 
for d = 2 and d = 3, respectively), A = 2tt^- 1 ^ 2 (1 - 
r 2 )a d - 1 /[dT(d/2)} (see e.g. [|), T(. . .) is the gamma- 
function, and d > 2 is the dimension of space, so that 
d = 2 corresponds to disks, and d — 3 to hard spheres. 
Furthermore, vq = {2u\/tt)~ 1 and kq = 4z/o in 2D, and 
v a = 5(3o- 2 V7F)~ 1 and k = 15z/ / 8 in 3D Equa- 
tions ([2j)- (J4j) differ from the hydrodynamic equations for a 
dilute gas of elastically colliding spheres only by the pres- 
ence of the inelastic loss term — ApT 3 / 2 which is propor- 
tional to the average energy loss per collision, ~ (1— r 2 )T, 
and to the collision rate, ~ pT 1 / 2 . 

It will be convenient for our purposes to rewrite 
Eqs. ©-((I]) in terms of the pressure p — pT, rather than 
the temperature. The energy equation Q becomes 



dp 
dt ' 
d 
dx 



dp 
dx 



-7P^-Ap 1 /2 p 3/2 



^0(7 





A set of hydrodynamic equations can be simplified if 
there is a time scale separation or, equivalently, a length 
scale separation, in the problem. For a freely cooling 
granular g basic time scale is the characteristic cool- 
ing time 



tc = 



(6) 



where po is the average gas density (the total gas mass 
divided by the volume of the channel), and po is a char- 
acteristic value of the initial pressure. There are two 
characteristic length scales related to t c . The first is the 
sound travel distance 



I, = 



7 



v/2 



Ap 



which is of the order of the distance a sound wave with 
speed c s = (7P0/P0) 1 / 2 travels during the time t c . The 
quantity L is the same as the length scale I introduced 
inRefs. il- 



The second characteristic length scale is the heat dif- 
fusion length 



1/2 
kqPq t 

3/2 



1/2 



.1/2 



A 1 / 2 



which, up to a numerical pre-factor, coincides with the 
critical length 




(7) 



predicted by the linear theory of the clustering instability. 
The ratio l s /l d is of order (koA)" 1 / 2 - (1 - r 2 )" 1 / 2 . As 
we have already assumed a strong inequality 1 — r 2 <C 1, 
this ratio is very large: l s /ld ^ 1- Throughout the rest 
of the paper we will also assume that the channel length 
L is much shorter than the sound travel distance l s . This 
hierarchy of length scales brings about a reduced set of 
equations, in much the same way as in hydrodynamics of 
optically thin gases and plasmas that cool by their own 
radiation [U, [H [H, Q . Note that the length scale 
separation L <g; l s is equivalent to a time scale separa- 
tion: the sound travel time through the channel, L/c s , is 
much shorter than the characteristic cooling time t c . As 
a result, sound waves rapidly make the pressure (almost) 
homogeneous throughout the channel. The subsequent 
slower evolution of the gas proceeds on the background of 
an almost homogeneous (but in general time-dependent) 
gas pressure, while typical Mach numbers of the flow are 
much less than unity. In a more formal language, this 
reduction of the hydrodynamic equations corresponds to 
elimination of acoustic modes. 

Before we perform the reduction procedure, let us in- 
troduce rescaled variables. We will measure the distance 
along the channel in the units of l cr , rescale time by t c , 
and measure the gas density, pressure and velocity in the 
units of po, po and l cr /t Cl respectively. Keeping the orig- 
inal notation for the rescaled variables, we observe that 
Eq. does not change, while Eqs. ([3]) and {5} become 





(8) 



dv 
dx 



IP^r- ~ 2p 1/2 p 3/2 



£2 (7 - 1 




where E\ = kqA/2, and £2 = VqA/2, and e\ ~ £2 ~ 
1 — r 2 -c 1. We will limit ourselves to the zeroth or- 
der approximation with respect to this small parameter 
and send £1 and £2 to zero. The continuity equation @ 
does not change. The momentum equation (jSJ becomes 
dp/dx — 0, therefore p = p(t) is independent of x. The 
energy equation becomes 



-7p(t) 



dv 
dx 



2p 1 ' 2 p(tfl 2 
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Pit) 



3/2. 



dx 



_L_ d_ 

yff) OX 



The rescaled length of the channel is 



( y5F/2)(l -r 2 ) 1 / a p g.L in 2d, 
y/lGn/75 (1 - r 2 ) 1 / 2 p cr 2 L in 3d. 



(10) 



(11) 



Note that, in the rescaled variables, the rescaled length 
of the channel £ coincides with the rescaled total mass 
of the gas, J p(x, t) dx. 

HI 



To get an explicit expression for p we integrate Eq. 
over the whole channel. Assuming either periodic, or no- 
flux boundary conditions (BCs) at the channel ends x = 
and x = C, we obtain 



Pit) 
p(t) 3 / 2 

where we have introduced the spatial average 

(...), = \[ (...)*. 



(12) 



For the low Mach number flow, Eq. (112| describes, in the 
leading order, the global energy balance of the gas, see 
Section [VTl C below. Equations ©, (JTUJ) and (JT2J for 
p(x,t), v(x,t) and p(t) make a complete set of reduced 
but fully nonlinear equations for the low Mach number 
flow of a freely cooling granular gas in a channel geom- 
etry. As is usually the case for low Mach number flows, 
the viscous terms dropped from the reduced formulation, 
while the heat diffusion term remains. 

The rescaled length/mass of the system £, see Eq. (TTTI) . 
is determined by the relative role of the inelastic energy 
loss and heat diffusion. As we will see shortly, C controls 
the main properties of the cooling dynamics. For com- 
parison, the characteristic initial pressure po only sets the 
time scale for the dynamics. To facilitate future compar- 
isons of the theory with MD simulations, we rewrite the 
parameter £ in a slightly different form: 



C 



^/l6TT{l-r 2 ) N p 



in 2d, 
in 3d . 



Here N p is the total number of particles in the channel, 
and L y and L z are the transverse channel dimensions. 



III. LAGRANGIAN DESCRIPTION AND 
NONLOCAL REACTION-DIFFUSION 
EQUATION 

Remarkably, it is possible to bring the three equa- 
tions (0) , (TIT))) and (fT2|) to a single nonlocal equation of a 
reaction-diffusion type. Let us first introduce Lagrangian 
mass coordinates [20|. It is convenient to choose a ref- 
erence frame so that v(x = Q,t) = 0. For the periodic 



boundary conditions (BCs) one can always achieve this 
by exploiting the Galilian invariance of the hydrodynamic 
equations to get rid of the center-of-mass motion. This 
sets v(x — 0, t) — 0, where x = is the center-of-mass 
coordinate. For the no- flux BC (impenetrable walls), a 
natural choice of x — is at one of the walls, where the 
gas velocity is again zero. Then a convenient choice of 
the Lagrangian mass coordinate is 



m(x,t) = / p(x ,t)dx , 



(13) 



which is simply the (rescaled) mass content between the 
Eulerian points and x. The inverse transformation is 



x(m, t) 



dm' 



(14) 



la P( m \ t) ' 

In the Lagrangian coordinates Eqs. ([2]) and (fT0| become 



8(1 

at \p 



dv 
9m' 

= - 1PP ^-2p^p^ 
om 

om \ om p 



(15) 



(16) 



As the total rescaled mass of the gas is equal to the 
rescaled channel length C, we define the spatial average 
in the Lagrangian coordinate as 



1 



and rewrite Eq. (TT2"]) as 

m 



P (t) 3 / 2 



= -2 



(. . .) dm . 



pV 2 (m, t) 



(17) 



It is convenient to introduce a new rescaled variable 
w(m, t) = p~ 1/>2 (m, t) and a new rescaled time 



T=- [ P l / 2 {t')dt' . 

7 Jo 



(18) 



Then, by eliminating d m v and p from Eqs. (|T5|) - (fT7| . we 
can reduce these equations to a single integro-differential 
equation of a reaction-diffusion type: 



d w 2 , \ 

w— — = —W + w (w) + 

OT 



dm 2 



(19) 



Equation (Tl9|) describes a broad class of slow Id flows in 
freely cooling nearly elastic granular gases. In particular, 
this equation encodes the development of the clustering 
instability: from a weakly perturbed HCS (after a brief 
acoustic transient) all the way to the strongly nonlinear 
stage. Indeed, let us rewrite Eq. (TiT]) in terms of the new 
variable w(m, r) and new time r: 



1 dp 
p(r) dr 



-27 (w(m, t)) . 



(20) 
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Once Eq. (jT9|) for w(m, r) is solved, we can calculate 
the pressure p(t) from Eq. (|2T))) and then return to the 
(rescaled) physical time t using Eq. (fl8|) : 



i = 7 



o P 1/2 (r') 



(21) 



Furthermore, using Eq. (|15p and the condition u(m = 
0, t) = 0, we can find the gas velocity: v(m, t) = 
J dtw 2 (m' ,t)dm! . Finally, we can return to the 
Eulerian coordinate by using Eq. (|14[) : x(m, t) — 
/ Q m w 2 (m' , t)dm'. 

Notably, equation (fill)) is parameter-free: the only pa- 
rameter entering the problem [except possible parame- 
ters introduced by the initial condition w(m, 0)] is the 
rescaled system length/mass £. Conservation of the total 
mass of the gas in the channel appears in the Lagrangian 
formulation as the conservation law 



(w 2 (m, t)) = 1 . 



(22) 



easily verifiable from Eq. (jT9 



IV. SIMPLE TESTS: HCS AND LINEAR 
THEORY OF CLUSTERING INSTABILITY 



As a first test of Eqs. (jl9|) and (I20|l . let us consider 
a HCS. Here at t — we have (in the physical units) 
pirn, t = 0) = pa = const, T(m, t — 0) = To = const 
and v(m,t — 0) = and, therefore, w(m,t = 0) = 1 
and p(t = 0) = po = poT Q . As the gas density remains 
constant in space at t > 0, we can rewrite Eq. (|19[) as 



dw(r) 
dr 



-1 + w 2 {t) 



(23) 



The solution of this equation with the initial condition 
w(0) = 1 is of course w(t) = 1: the gas remains spatially 
homogeneous. Now we use Eq. (fl~7|) and obtain 



P(t) 
p{tf/ 2 



(24) 



which yields, in the rescaled variables, Haff's law for the 
gas pressure: 



P(t) = 



1 



(25) 



The next test of Eq. (119|) is the linear stability anal- 
ysis of a HCS. While the reduced Eq. (|19p is not sup- 
posed to capture the evolution of small perturbations 
with an arbitrary polarization, it must reproduce cor- 
rectly the evolution of the clustering mode in the limit 
when the perturbation wavelengths are small compared 
with the sound travel distance l s . Let us show it to be 
indeed the case. We look for the solution of rescaled 
Eq. (|19p in the form w(m,r) = 1 + Sw(m,T), where 
\5w(m, t)\ <C 1. [Correspondingly, the rescaled density 



perturbation Sp(m, r) = — 2Sw(m, r).] One can represent 
5w(m, t) as a linear superposition of sines and cosines of 
km with different (rescaled) wave numbers k. This fact, 
in conjunction with the BCs at the ends of the channel, 
guarantees that (5w(m,T)) = 0. Then Eq. (fT9]l yields 

d d 2 

— Sw(m, t) — Sw(m, r) + — — ^Sw(m, r) . (26) 
or dm 1 

For a single mode perturbation with wave number k we 
obtain 



(27) 



(28) 



Sw(m, t) = 5w(m, 0) e TkT 
with the growth/damping rate 

f fc = 1 - k 2 . 



For k < k* = 1 (correspondingly, k > fe* = 1) Eqs. (|27|) 
and (|28|) describe an exponential growth (correspond- 
ingly, decay) of a small single-mode perturbation in time 
r. Recalling that we rescaled the coordinate to the crit- 
ical length l cr , provided by the complete (unreduced) 
linear theory, we immediately notice that Eq. (f2"8")) cor- 
rectly predicts the instability threshold. To go back to 
the physical time t we substitute, in the leading order, 
Haff's law (|25[) into Eq. (jTSj) and obtain, after elementary 
integration, 



t = - In (1 + t) 

7 



(29) 



Plugging it into Eq. (|27|) . we obtain an algebraic growth 
of the small perturbations in the physical time: 



5w{m, t) — Sw(m, 0) (1 + t) 



r fc /7 



(30) 



The growth exponent T = Tk/j, with Tk from Eq. (|28|) 
coincides with that obtained from the complete linear 
stability analysis Q , if we assume there kl s ^> 1 (in the 
physical units) and consider the clustering mode, rather 
than the two decaying acoustic modes. Figure [1] shows 
this comparison in a graphic form. At kl s < 1 the iso- 
baric growth rate underestimates the true growth rate, 
but in the region of kl s » 1 excellent agreement is ob- 
served. The comparison with the complete linear stabil- 
ity analysis is instructive for two more reasons. First, as 
was observed by McNamara 0], for kl s ^> 1 the pres- 
sure perturbations of the clustering mode vanish in the 
leading order in l/(kl s ). That is, the linear density and 
temperature perturbations grow on the background of 
an (almost) constant pressure. Second, the viscosity ef- 
fects do not affect the growth exponent in this regime 0] • 
As our reduced formalism shows, the last two properties 
persist, for the low Mach number flow, in the nonlinear 
regime as well. 

Having successfully tested our reduced model in these 
two simple cases, we now consider nonlinear evolution. 
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FIG. 1: (Color online.) The growth exponent Tk of the clus- 
tering mode versus the rescaled wave number kl s (a) and ver- 
sus log (kls) (b) as predicted from the complete linear stability 
analysis of a HCS [f| (the thick black line) and from the re- 
duced equation (|19p (the thin red line). The physical (not 
rescaled) units are used, and the parameters are 7 = 2 and 
k cr l 3 = 100. At kl s ^ 1 the isobaric growth rate ()19[) under- 
estimates the actual growth rate, but in the region of kl s ^> 1 
excellent agreement between the two results is observed. 



V. NEGLECTING HEAT DIFFUSION CAUSES 
A DENSITY BLOWUP 

As we mentioned earlier, the only governing parameter 
in Eq. (|19p , except parameters introduced by the initial 
condition, is the rescaled system length/mass £. In the 
limit of £ 3> 1, and for a sufficiently large-scale initial 
condition, one can drop the diffusion term in Eq. (|19[) . 
This approximation is valid as long as the solution re- 
mains large-scale. At the level of linear stability analy- 
sis this (intermediate-wavelength) approximation is fully 
justified. Here Eq. (|2"5|) becomes Tj~ — 1, and one is 
interested in the nonlinear development of the growing 
perturbations. With the diffusion term dropped we ob- 
tain 



dw 



-1 + w (w) 



(31) 



This nonlinear integro-differential evolution equation is 
exactly soluble for any initial data w(m, 0). The com- 
plete solution is presented below. The main result here 
is that, for any inhomogeneous initial condition, the so- 
lution of Eq. (|3Tj) develops a zero w (hence an infinite 
density) in a finite time. Let us first discuss the proper- 
ties of the solution in a close vicinity of the singularity 
w — ► 0. In the leading order we can neglect the inte- 
gral term in Eq. (f3Tj) and obtain dw/dr — —1, so that 
w(m,T) — w(m) — t, where w(m) is a smooth function. 
The singularity occurs in the Lagrangian point Too that 
corresponds to the minimum of w(m). The leading order 



behavior of the (rescaled) gas density near the singularity 
is described by the following equation: 



P(to,t) = 



Tr — T + - - 



2 dm' 



(to ) (to - mo) 



(32) 



where the time of singularity r c = w(m ). The singu- 
larity structure, as described by Eq. (j32|) . coincides with 
that exhibited by a family of exact solutions of the full 
set of ideal hydrodynamic equations [that is, Eqs. ([2])- 
((!]) without the viscous and heat diffusion terms], re- 
ported in Ref. [ij], Hq |- At r = t c the density blows 
up as p(to,t c ) ~ (to — mo) -4 . Going back to the Eule- 
rian coordinate, we obtain a finite-mass density blowup 
p(x,T c ) ~ \x — a; |~ 4 / 5 , where Xq is the Eulerian coordi- 
nate of the singularity. We refer the reader to Ref. [HI 
for a detailed analysis of the structure of this singularity, 
as observed in the gas density, temperature and velocity. 
Notably, the pressure field does not have any singularity 
in the exact solutions [l7l [l8| , and is approximately con- 
stant in a narrow region around the density singularity. 
That is, the density blowup, as featured by the exact so- 
lutions of ideal granular hydrodynamics [Tt], EH , locally 
obeys an isobaric scenario, as was noticed in Ref. [18l | . 
This provides the reason why the same type of singularity 
appears in our reduced low Mach number theory. 

Now we present a complete solution of Eq. (f3"Tj) . First, 
we obtain a closed evolution equation for the (necessarily 
positive) quantity x( T ) — (w(m, r)) by integrating the 
both sides of Eq. (|3"Tj) over to from to C: 



dx(r) 
dr 



(33) 



We consider the solution of this equation with the initial 
condition 



Xo = (w(m,0)) = (p(m,0) 



-1/2 



< 1 



The solution can be written as 

Xo - tanh(r) 



X(r) 



1 - Xo tanh(r) 
Now we can rewrite Eq. (f3"Tj) as 

— - x(T)w(m,T) = -1 , 

where x( T ) is given by Eq. ([M)) . Equation (f?5|) is easily 
soluble: 



(34) 



(35) 



u>(to, t) 



w(m, 0) + xo [cosh (r) — 1] — sinh (r) 
cosh (r) — xo sm h (r) 



(36) 



The presence of the factor xo [cosh (t) — 1] — sinh (r) in 
the numerator of Eq. ()36|) causes, for any (non-constant) 
initial data w(m, 0), a singularity w — > in a finite time. 
The singularity occurs at the Lagrangian point mo where 
the function w(m, 0) has its minimum, at time 



In 



Aw + (Aw 2 + 1 - Xo) 
1 - Xo 



1/2 



(37) 



7 



where Aw = w(mo,0) — X o- Note that Aw = 
(l/£)/ [w(mo, 0) — w(m, T)]dm < 0. 

Now we compute the (rescaled) pressure p(r) from 
Eq. (T20|) [note that the right hand side is simply x( T ) 
given by Eq. l[34|)]. 



p(r) = (cosh r — xo sinh r) 



2- 



(38) 

and use this result in Eq. (|2 1 1) for the rescaled physical 
time: 

dr' 



i = 7 



(coshr' — xo sinnr') 



(39) 



For 7 = 2 (a 2D gas of disks) this integral is elementary 
and the result is 



t 



cothr - xo 
Now we can express r through t, 



t = arccoth [ — + Xo 



(40) 



(41) 



and rewrite Eqs. (J2U) and (f3"8")) (for 7 = 2) as 

w(m,t) = ~{[w(m, 0) - xo] ^4 + 4 Xo t - (1 - x§)i 



+ 2xo-(l-Xo)*}- 



and 



P = 



1 + Xo* - (1 - Xo) T 



(42) 



(43) 



So, the solution for 7 = 2 is surprisingly simple. We 
remind that, in view of the chosen rescaling, the initial 
condition w(m, 0) must obey (u> 2 (m, 0)) = 1. To return 
to the HCS and Haff's law in Eqs. gU) and J43} one 
should put there w(m,0) = xo = 1- Equation (|43[) shows 
that Haff's law is an upper bound for the thermal energy 
loss rate: any deviation from homogeneity brings about 
Xo < 1 an d a slower thermal energy decay. 

Let us note that the solution (|3"4"} for x( T ) vanishes 
at t* = (1/2) ln[(l + Xo)/(l — Xo)] an d becomes nega- 
tive at larger r. This is in apparent contradiction with 
the positivity of w that dictates x( T ) ~ (w(m,T)) > 0. 
The contradiction is resolved by noting that r* is always 
greater than the singularity time r c , beyond which the so- 
lution does not apply. [To see that r c < t* one can use, 

inEq. that Aw+(Aw 2 + l- X 2 o) 1/2 < (l-xlf' 2 
for any Aw < 0.] Similarly the pressure as predicted by 
Eq. (|38|) or Eq. (|43|) would start increasing at some time. 
At physically meaningful times t < r c , however, we have 
x(t) > 0, and the pressure always decreases in accord 
with Eq. flU}. 

As a simple illustration of our solution (f36|) . let 
us chose the following initial condition: w(m, 0) = 
[1 + 5cos(2nm/C)} 1/2 , < S < 1. In this case 



1 

X = - 

7T 



vT^E 
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FIG. 2: The density history of a freely cooling gas of inelas- 
tic hard disks in a 2D channel in the zero-heat-diffusion limit. 
The rescaled initial density p(m, 0) = [1 + 0.1 cos(2-7rm/£)] -1 . 
The rescaled system length/mass C = 100. Panel a: the 
rescaled inverse density of the gas, 1/p, versus the Lagrangian 
mass coordinate m at times r = 0, 1.5, 2.5 and the time of 
singularity r c ~ 2.8755. Panel b: the rescaled Eulerian coor- 
dinate x versus m at the same times. Panel c: 1/p versus x 
at the same times. The sequence of curves is self-explanatory. 



where E(. . .) is the complete elliptic integral of the sec- 
ond kind, see e.g. [30| ■ Figure [2^ shows, at different 
times, the rescaled inverse density l/p(m, r), as obtained 
from Eq. $55} , for 6 = 0.1 and C = 100. Figure [2Jd de- 
picts, at the same times, the rescaled Eulerian coordinate 
x = f™w 2 (m',T)dm' versus the Lagrangian coordinate 
to. Figure shows the rescaled inverse density in the 
rescaled Eulerian coordinates and illustrates the emer- 
gence of the cusp density singularity at x — C/2. The 
inverse density behaves like (m-mo) 4 at small to — too in 
the Lagrangian coordinate, and like (x — xq) 4 ^ 5 at small 
x — xq in the Eulerian coordinate. This simple example 
is instructive as, for <5 <C 1, this initial condition corre- 
sponds to a small single-mode density perturbation, so 
the initial evolution is describable by the linear theory. 
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VI. HEAT DIFFUSION ARRESTS THE 
DENSITY BLOWUP 

A central result of this work is in that, no matter how 
small initially, the heat diffusion term in Eq. (Til?)) arrests 
the density blowup. An emerging balance of the inelas- 
tic cooling and heat diffusion leads to existence of steady 
state solutions of Eq. p^|) . These solutions describe novel 
cooling states of the granular gas, where the (inhomoge- 
neous) density profile is time-independent, while the (ho- 
mogeneous) pressure continues to decay with time. We 
found that, in our rescaled variables, the density profile 
of the novel cooling state is uniquely defined by the pa- 
rameter £,. For sufficiently large values of the rescaled 
length/mass, C ^> 1, the maximum gas density of the 
novel cooling state is exponentially large in C In the 
low Mach number theory, considered in this work, the 
novel cooling states represent global attractors, as they 
develop for any inhomogeneous initial conditions. Fi- 
nally, the novel cooling states represent exact solutions 
of the complete, unreduced set of hydrodynamic equa- 
tions ©-©• 



that, in virtue of Eq. (|46|) , is obeyed automatically for 
the periodic or no-flux BCs. 

Equation (|46|) has appeared in numerous applications, 
and its solutions are well known. It is convenient to in- 
terpret / as a coordinate of a Newtonian particle of unit 
mass, moving in a potential U(f) — / 3 /3 — / 2 /2. The 
"total energy" E is conserved: 



E= I (df_Y , f P 
2 [dm. 



(49) 



For the bounded (spatially oscillating) solutions, — 1/6 < 
E < 0, we can write 



/ 3 P 

- E 

3 2 



(f-a)(f-b)(f-c) 



(50) 



where a > b > c are the real roots of the cubic poly- 
nomial. Then the bounded solutions of Eq. ([46]) can be 
written as 



/(to) = c + (a — c) dn 2 




(51) 



A. Steady state density profiles 

Steady-state solutions of Eq. (TIT)]) are described by the 
equation 



where 



d 2 w 
dm 2 



w — (w) w 2 



(44) 



Notice that, although obtained from our reduced, low 
Mach number theory, Eq. (TJ3|) also follows from the full 
set of hydrodynamic equations ©-(Q, if one assumes a 
homogeneous pressure and zero fluid velocity, and trans- 
forms to the Lagrangian coordinates. 

Equation (|4"4")) is defined on the interval < to < C, 
at the ends of which we demand either periodic, or no- 
flux (zero first derivative) BCs. The solutions we are 
interested in must obey the conservation law ([2"!?]) . To get 
rid of the (a priori unknown) factor (w), we introduce a 
new variable 



/(to) = (w) w(m) 



and obtain 



(45) 



(46) 



(47) 



The conservation law (|22|) enforces a normalization con- 
dition 



#4-/+/ 2 =o- 

dm 

Once / is found, one can restore w via 

/ 

w = 



</ 2 > = (/> 



(48) 



(52) 



and dn is one of the Jacobi elliptic functions, see e.g. [301 ] . 
There are two limits when Eq. (|5Tj) simplifies. In the limit 
of E = -1/6 + 6E, < &E< 1, the solution, /(to) = 1 + 
y/2SE cos to, corresponds to a small-amplitude sinusoidal 
modulation of the HCS w(m) = 1. In the limit of E — ► 0, 
we have a = 3/2 and b = c = 0, so that 



/(TO,£^0) = ^dn 2 (^,l) 



= - cosh 



(?)■ < 53 > 



Using Eqs. (j47|) and (j5lj) , we rewrite the steady state 
solutions in terms of w(m): 



w(m) 



c + (a — c) dn 2 (y^p m , s) 



c+(a-c)|j 



(54) 



where K(s) is the complete elliptic integral of the first 
kind. The lagrangian spatial period, or wavelength, of 
the solution (|54|) is 



n 



24 



K( S ). 



(55) 



In the limit of s — > (or E — > —1/6), the wavelength |55|) 
reaches its minimum value 2tt. If the rescaled channel 
length C is less than 27T (for the periodic BCs), or less 
than 7r (for the no-flux BCs), the only possible steady 
state is the constant density state w{m) — 1 correspond- 
ing to Haff's law. This result is in full agreement with the 
linear stability analysis of Eq. (Til)]) , see Eq. (|28p . When 
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C exceeds 2ir (for the periodic BCs), or ir (for the no-flux 
BCs), the HCS bifurcates into an inhomogeneous steady 
state (|54|) . In general, the rescaled channel length/mass 
£ must be equal, by virtue of the BCs, to an integer num- 
ber of II (for the periodic BCs), or to an integer number 
of 11/2 (for the no-flux BCs). For sufficiently large value 
of C, therefore, a whole family of steady state density 
profiles exists. Which of the steady state solutions is 
selected by the dynamics of Eq. (fl9|) ? 



the 



Selected steady-state solutions: 
inhomogeneous cooling states 



We performed extensive numerical simulations with 
Eq. (|19[) . using a specially developed numerical scheme 
described in Appendix A. Both periodic, and no-flux BCs 
were used. We observed that, when < L < 2ir (for the 
periodic BCs), or < C < ir (for the no- flux conditions), 
the HCS appears, as expected. When C exceeds 2n (for 
the periodic BCs), a weakly inhomogeneous steady state 
density profile sets in. As C increases further, the weakly 
inhomogeneous states develops into a strongly inhomoge- 
neous states. The simulations showed that the rescaled 
length/mass of the gas, C, uniquely selects the emerging 
steady state density profile, while the initial w-profile 
does not play any role in the selection. For a given C 
the dynamics always selects, out of the family of steady 
state solutions (|54[) . the one with the maximum possible 
wavelength II: 



C = 



II for the periodic BCs , 
II/2 for the no-flux BCs . 



(56) 



Snapshots from a typical simulation (one of many that 
we performed) for the periodic BCs are shown in Fig. [3] 
The initial condition is this example was 

w 2 (m,0) = 1 - 0.1 cos(27rm/£) - 0.15 sm(2irm/£) 
+ 0.2 cos(47rm/£)- 0.05 sm(47rm/£). (57) 

The rescaled system length/mass C — 50 was sufficiently 
large to fit in steady state solutions with several oscil- 
lations. Nevertheless, the dynamics selected the solu- 
tion with the spatial period equal to the rescaled system 
length C. 

Figures [4] - [6] depict our analytical solutions (|54|) in 
the Lagrangian coordinate, and the corresponding den- 
sity profiles in the Eulerian coordinates, for three dif- 
ferent values of the parameter C Here we assumed the 
periodic BCs and (arbitrarily) chose the position of the 
minimum of w{m) to be in the middle of the channel. 

The maximum (rescaled) gas density versus the 
rescaled channel length C, predicted by Eqs. ([54")) and 
(|55|). is shown in Fig. [7l This dependence can serve as 
a bifurcation diagram of the system. One observes, at 
C > 27r, a supercritical bifurcation from the HCS to an 
ICS. 
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FIG. 3: Numerical ui-profiles at times r = 0, 2, 4 and 72, and 
1/w-profile at time r = 72, for C — 50 when starting from the 
initial condition (|57|l . The two panels for r = 72 also show, by 
circles, the single hole asymptotes (|60[1 and (|6ip . respectively. 
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FIG. 4: The inhomogeneous cooling state for C = 7.025. 
Panel a: the Lagrangian steady state solution w(m) as pre- 
dicted by Eq. (|54[1 . Panel b: the rescaled steady state gas 
density p versus the rescaled Eulerian coordinate x. 



One can see that, as the parameter C increases, the 
maximum gas density in the cluster grows very fast [note 
that Fig. [6]b shows the density in logarithmic scale] . Let 
us consider the asymptotic form of the solution at C 3> 1 
in some detail. The density maximum in this case is 
exponentially large (3lj . This is due to the behavior of 
the s — ► 1 asymptotics of the steady-state solution, see 
Eq. In this case the "energy" E is very small, and 

can be expressed through the rescaled system length as 
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FIG. 5: Same as in Fig. H but for C = 8.886. 
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FIG. 6: Same as in Fig. H but for C = 19.869. Notice the 
logarithmic scale in panel b. 



\E\ ~ 72 exp(— C). The maximum value of w(m) is 




(58) 



To obtain the minimum value of w(m) (that corresponds 
to the maximum value of the density) , it is convenient to 
use the exact relation w m i n — b/ yJJJ) and calculate the 
asymptotic value of b at \E\ -C 1, or L 1. The result 
is 

w min - V24Ce~ c/2 , (59) 
By virtue of Eq. (|53p , the asymptotics of the steady state 
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FIG. 7: The bifurcation diagram of the freely cooling granular 
gas in a channel. Shown is the maximum (rescaled) steady 
state density of the gas versus the rescaled channel length C, 
predicted by Eqs. (|54[) and (|55[) . Panel b focuses on a vicinity 
of the supercritical bifurcation point C = 2tt. 



solution ([54]) at \m\ < C/2 is 



w (m) 



cosh 2 



(?)■ 



(60) 



where, for convenience, we have written the solution on 
the interval —C/2 < m < C/2 and used the approxi- 
mate equality (/) ~ Q/C. To calculate the asymptotics 
of Eq. ([51)1 at \m\ > 1, we can deal directly with Eq. [(IB"]) 
and neglect the / 2 term. The solution of the resulting 
elementary equation is a linear combination of e m and 
e _m . The two arbitrary constants can be determined 
from the two conditions at \m\ = C/2: df /dm = and 
w o = f I VTT) = w min, where w mm is given by Eq. ((59]). 
We obtain 

w (m) ~ V24Ce" £/2 cosh(£/2 - |m|) , |m| > 1 . (61) 

Note that the asymptotes (pO")) and (|6T|) coincide in their 
common region 1 <C |m| <C >C/2, where each of them 
yields 



w (m) 



(62) 



Note that (wo(m)) ~ y/6/C is determined by the asymp- 
tote (|60|) . We compared the asymptotes (|60|) and (|6"T]) 
with the numerical solution, shown in Fig. [3l at a late 
time r = 72. Employing the periodic BCs, we shifted 
the numerical solution in m so that the maximum of 
w(m, t — 72) is at m = 0. One can see that the agree- 
ment is excellent. 

As higher w corresponds to a lower gas density, the 
region of the maximum of w corresponds to a hole in the 
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density. Therefore, we will call the approximate solution, 
fully determined by Eqs. ([SO]) and (fBTj) . the hole solution. 
The rescaled steady-state gas density, in the limit of £ 
1, is 



This form is convenient in the vicinity of m = £/2. The 
case of to = —C/2 can be treated similarly, and the ex- 
pressions that follow are valid in both cases. For \m\ ^> 1 
Eqs. dHHJ) and §T§ yield 



p{m) 



3£ 



cosh 



— J , |m| < C/2 , (63) ~ — — 6 £e~ £ [£ - 2|to| + sinh(£ - 2|m|)] . (69) 



and 



Pi m ) - cosh 2 I - - \m 



C 



Iml » 1 



(64) 



and the maximum and minimum density values are 



Prnax — ' Prnin — ^£ 



(65) 



Note that Eqs. (|58|) - ([65)) work very well already for mod- 
erate values of £. For the dilute hydrodynamics to be 
still valid in the gas density peak region, we must de- 
mand that the peak density be much less than the close 
packing density. In view of the exponential growth of the 
maximum density with the parameter £, see Eq. (f6"5"|) . 
this leads to a stringent condition: 



po<r d < 24£ e 

If this condition is not fulfilled, the dilute theory will 
break down, and the attempted density blowup will be 
regularized by close-packing effects. 

The general form of the steady state density profile in 
the Eulerian coordinates is quite cumbersome. However, 
its asymptotic form at £ 1 that corresponds to the 
Lagrangian profiles (f6"0"|) and (fBTj) is both elementary and 
instructive. For Eq. (fBT))) one finds, after some algebra, 



( ^ / 3£ 
wq{x) = W — cos 



1 / 8x 2 

-arccosfl- — 



This asymptotics is valid at e~ c <C 1 — 2\x\/£, that is al- 
most over the whole channel \x\ < C/2 except in a narrow 
region. This region, however, includes a significant part 
of the gas mass, as evidenced by the size of this region 
in the Lagrangian coordinate and by the non-integrable 
diverging power-law asymptotics of the gas density: 



p{x) 



1 



£-21. 



at 



e- c < 1 



2N 

c 



< 1. 



(67) 



There is of course no actual density divergence here, 
as Eq. (|67|) does not hold close to the end points: at 
1 — 2|ir|/£ < e~ c . To find the density profile in this 
exponentially narrow region, we express the relation be- 
tween x and to as 



w (to') dm' 

C/2 rC/2 

WQ(m')dm— / w 2 (m')dm 



rC/2 

= C/2- I wl(m')drri . 



(68) 



Equations (|64[) and (|69p determine, in a parametric form 
and in elementary functions, the density profile in the 
region sufficiently far from the density minimum. Still 
simpler results can be obtained in the following two sub- 
regions. The first is the common region C/2 — \m\ ^> 1 
but \m\ > 1. The asymptotics of Eqs. ((Ml) and (JM]) 
at £/2 - \m\ > 1 become p = e 2 l m l/(6£), and \x\ ~ 
£/2 - 3£e- 2 l™l, therefore p = (£ - 2|a;|)- 1 which coin- 
cides with the asymptotics (|57)) of Eq. (|6^|) . The second 
limit corresponds to C/2 — \m\ <§; 1. Here Eq. (|64|) be- 
comes 



Fy ' 24£ 



whereas Eq. ^ yields \x\ = C/2 - 24£e~ £ (£/2 
The resulting Eulerian density profile is 



i-C— )'(-- 

\24CJ \2 



H). 



(70) 



C. Energy decay for the ICSs 

Now let us consider the evolution of the (rescaled) total 
energy of the gas, 



E tot (t) 



P 



pv 



7-1 



dx . 



(71) 



where the first term under the integral is the thermal 
energy density, and the second term is the macroscopic 
kinetic energy density. For the low Mach number flow, 
that we are dealing with in this work, the first term is 
almost independent of x, while the second term is neg- 
ligible. As a result, the energy decays, in the leading 
order, in the same way as the pressure. The pressure de- 
cay p(r) is described by Eq. (j2TI)) . whereas to go back to 
the physical time we use Eq. (f21~j) . For our steady state 
solutions we arrive at a generalized Haff 's law 



Pit) 



1 



(i + (w)ty 



(72) 



As (w) < {w 2 ) 1 / 2 = 1, the energy decay for the ICS is 



always slower than for the HCS, see Eq. ([25j) . A more 
explicit form of the generalized Haff 's law (f72l) is 



Pit) 



1 + tt c+(a-c) 



K( S ) 



(73) 



12 



Now we consider the particular case of the single hole 
solution Wo(m). As (wo(m)) ~ y/6/C, we obtain for the 
pressure (in the physical units) 



(t) = p Q cxp (-2VQiC~ 1/2 t 



(74) 



Using Eq. (f2Tj) . we find the original (physical) time t in 
terms of r (again, in the physical units): 



This yields a generalized Haff's law 

Pa 



Pit) = , 

(1 + t/t c 

with a characteristic cooling time 



(75) 



(76) 



(77) 



As £ 3> 1, the cooling time f c is much longer than the 
cooling time t c corresponding to the HCS: 



1/2 



< 1. 



(78) 



D. Relaxation to the single hole state 

Here we study the late-time dynamics of relaxation of 
the cooling gas towards the single hole state: the cooling 
state observed for C 3> 1, that is, for l cr -C L <C l s . We 
put w(m,T) — wo(m) + wi(m,r), where wo(m) is the 
single hole asymptotics (J5DJ), and linearize Eq. (TT9"]) with 
respect to the small correction w\. We obtain 



dwi , . . . . . 2 <9 2 u>i 

w o-g^r = ( 2w o (wo) - l)u>i + (toi) w + 



dm 2 



(79) 



In the language of the linear stability analysis, the con- 
servation law (|22p becomes (w;o(to)u;i(to, t)) = 0. Inte- 
grating Eq. (|79p over the box, one can see that, once this 
condition holds at r = 0, it continues to hold at r > 0. 

As will become clear shortly, a natural complete set of 
eigenfunctions for the linear equation (|79|) is provided by 
the following eigenvalue problem: 



2/«( m ) + [-1 + Kw Q {m)} y n {m) = . 



(80) 



for the eigenfunctions y n (m) obeying the BCs y n (±oo) = 
0. (Here we have moved the boundaries to infinity which 
is accurate with an exponential accuracy in the large pa- 
rameter C ^> 1.) Equation (|80|) can be viewed as a sta- 
tionary Shrodinger equation (with h = 1) for a particle 
with mass 1/2 and a fixed energy —1 in the Poschl- Teller 
potential well, see e.g. Ref. j33|. The depth of the well is 



determined by the eigenvalues A„ . The spectrum of this 
problem is discrete: 



(n + 2)(n + 3) 



71 = 0,1,2,3,. 



(6£)V3 

For even values of n one obtains even eigenfunctions: 
yT n {m) = A„cosh 2n+3 (|) 



(81) 



X 2*1 



5 1 

2'2' 



sinh' 



(?) 



(.82) 



whereas for odd values of n one obtains odd eigenfunc- 
tions: 

y?\m) = S rl cosh 2 "+ 4 (!)sinh(!) 



x 2 Fi 



A 7 3 . o f m \ 
n A — ,n A — ; — ; — sinh — 
2 2 2 V 2 / 



(83) 



Here 2-F1 (• ■ ■) is the hypergeometric function, and A n and 
B n are constants that we fix using the orthonormality 
conditions 



+00 



yk(rn)y n (m)w {m) dm = 6 kn , 



(84) 



the Kroneker delta. The fundamental mode yo(m) is 
even, it is proportional to wo(m): 

y a (m) = C w (m) = (75/2£) 1/4 cosIT 2 (to/2) , 

where Co = 5 1 / 2 6 _1 / 4 £ -3 / 4 . The next mode is the first 
odd eigenfunction yi(m), proportional to dwo(m)/dm: 



yi (to) 



Z 1 ^^ cosh -2 (to/2) tanh (to/2) 



2 7/4£l/4 

The next one is the second even eigenfunction 

3 3 / 4 Vb (3 cosh to - 4) cos!T 4 (to/2) 



2/2 (to) 



2 7/4£l/4 



and so on. Let us expand u>i(to,t) in this complete set 
of eigenfunctions: 



Wi (to, r) = ^ a n (t) y n (to) . 



n=0 

substitute this expansion in Eq. (|79p . multiply the re- 
sulting equation by yk(m), k — 0, 1, 2, . . . and integrate 
over to from —00 to 00. Using Eq. (|80[) . we arrive at the 
following equations for the time-dependent amplitudes 
Ofc(r): 



rfafc(r) 

dT 



-T k a k {r) for k ^0, 



and 



da (T) 
dr 



2 (w ) ao(r) + — V a 2n (r) (1/2,, 



(85) 



(86) 
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Here 



Tfc = A fe -2 (w ) 



(fc- l)(fc + 6) 
(6£)V2 



1,2,. 



(87) 



and we have used the equality Ao = (wq). The amplitude 
equations ([85]) and (|86)). together with the initial condi- 
tions afc(O), k = 0, 1, 2, . . ., enable us to solve the initial 
value problem for the evolution of the small perturbation 
wi(m, t). Equations (|85|) show that each of the odd and 
even modes k — 1,2,3,... evolve independently of other 
modes: the k = 1 mode has a zero decay rate (which is 
expected, as it is a translational mode) , while the higher 
modes decay exponentially in time r : 

flfc(r) =a fe (0)exp(-r fc r), k = 1,2,3,.... (88) 

The k — mode behaves quite differently from other 
modes, as it is affected by the rest of the even modes of 
the system, see Eq. (|86|). The solution of Eq. (j86[) is: 



^ OO 



O2n(0) (U2r, 



exp (2 (w ) t) 



1 y> «2n(0) (j/2n) 

Go ^ 



A 2 



A2n 

exp (-r 2n r) . (89) 



Now we prove that the term in the square brackets van- 
ishes. At r = the conservation law ((22)) can be written 
as 

( woirn) ^2 a 2n(0)y2n(m) j = 0, 



n=0 



which yields 



a (0) 



^ OO 

— ^2 a 2n (0) (w Q y 2n ) = . 



Co 



(90) 



By virtue of the identity (?/2n) = A 2n («W2n) [which read- 
ily follows from Eq. ([50])] . the left side of Eq. (J9Q1) coin- 
cides with the term in the square brackets in Eq. ([59")) . 
Therefore, the final result for ao(r) is 



^ OO 



a2n(0) (y 2n ) 



n=l 



A2n 



exp (-r 2n r) 



(91) 



ao(r) can behave non-monotonically at short times. 
However, it always decays at long times, and the domi- 
nant decay rate, atr> C 1 /' 2 , is T 2 - 

Figures [8] and [9] present a comparison of the linear sta- 
bility analysis with the simulation shown in Fig. [3J Fig- 
ure [8] shows, at late times, the deviation of the numeri- 
cal solution from the theoretical single-hole steady state 
asymptotics (|oT))) for the simulation shown in Fig. [3] As 
time proceeds, the deviation tends to zero as expected. 
Figure [9] compares the numerically observed decay rate 
of the deviation with the analytical result (|87[) for the de- 
cay rate T 2 that dominates at late times, and very good 
agreement is observed. 
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FIG. 8: The difference wi(m, r) between the time-dependent 
solution and the single-hole steady state asymptotics (|60l l at 
different (late) times for the simulation shown in Fig. [3] 



Using Eqs. ([75]) and ([57]) . we can see that the expo- 
nential decay in r of each of the eigenmodes k — 1,2,..., 
see Eq. (|55)l . becomes a power-law decay in the physical 
time: 



a k (t) = a k (0) [l + j 



(fc-l)(fc + 6) 
6-/ 



1,2,. 



with t c from Eq. (j77|) . The zero mode dynamics ([91]) can 
be represented as a superposition of terms, each of which 
decaying as a power law in the physical time. There- 
fore, the mismatch w(m,t) — wo(m) between the time- 
dependent solution w(m,t) and the single hole solution 
wo(m) decays, at long times, as ~ (t/t c )^ 4 ^ 3l \ 

Before concluding this section we note that the k — 1 
mode turned out to be marginally stable because we ne- 
glected corrections exponentially small with respect to 
the rescaled system length L. In a more accurate treat- 
ment this mode would cease to be a translational mode 
and acquire a non-zero (although exponentially small) 
damping rate in time r. This would lead to a power law 
decay of this mode in time t with a power exponent that 
is exponentially small in C 
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T 

FIG. 9: Testing the linear stability analysis of the single 
hole solution. The circles show, in the logarithmic scale, 
max m \wi(m, r)| (see Fig. EJ versus time r. The solid line 
depicts our theoretical prediction for long times, when the re- 
laxation is dominated by the mode 1/2, so that max m |u>i| = 
c exp(-r 2 T), where T 2 = 8/V(i£ ~ 0.462. The adjustable 
parameter Co = 33.5. 
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FIG. 10: Nucleation and coarsening of holes when starting 
from a small amplitude "white noise" density perturbation 
around w = 1. Shown is a small fragment of the system of 
rescaled length/mass C = 10 6 at indicated times. 



VII. COARSENING DYNAMICS AND 
STATISTICS OF HOLES 

Numerical simulations with Eq. (|19p show that, for 
a sufficiently large rescaled length/mass of the system, 
C >■ 1, many peaks of w (hence, holes of the gas den- 
sity) nucleate in the system [3l}. The nucleation stage, 
as observed numerically, is shown in the upper left panel 
of Fig. [TQJ The initial condition w(m, r = 0) simulated 
white noise, as we chose w 2 (m,T — 0) to be equal to 1 
plus a sum of a very large number of Fourier harmonics 
with (very small) random amplitudes drawn from a uni- 
form distribution. As evidenced by Fig. QJ2 the further 
evolution of the holes resembles Ostwald ripening [34| . 
At this stage nucleation of new holes does not occur any- 
more, and a competition between the holes begins. Un- 
derdense holes release their material into the environment 
and become more pronounced (even less dense), while 
holes with more material continue to suck the material 
in until they disappear. At some stage the holes which 
gas density previously decreased, reverse the trend and 
begin to densify. At the end of this coarsening process 
only one hole (that was the least dense in the beginning) 
remains and forms the single-hole solution (|60|) and (|6Tj) 
|32| . Clearly, the holes compete non- locally: via the spa- 
tial averaging term of Eq. Q19p . 

Can one build upon the analogy with Ostwald ripening 
and develop an asymptotic theory of the hole coarsening 
dynamics? Consider a late stage of the dynamics when 
there are N holes, located sufficiently far from each other, 
and centered at points mj, i = 1, 2, . . . , N. A simple the- 
ory assumes that the spatial shape of each hole coincides 
with that of the limiting steady state asymptotics (f60|) . 
but with its own amplitude Ai(t) that depends on time. 
The latter assumption is based on a remarkable fact that, 
up to exponentially small corrections, Eq. (|19p admits the 



following ansatz: 

w(m, r) = £ Mr) cosh" 2 ^-^ . (92) 
i=i ^ ' 

Plugging it into Eq. (|19|) and neglecting exponentially 
small overlap terms, we find that the equation is satisfied 
once the following N relations hold: 

Mr) = S(t)Mt) ~ \ , t = l,2,...,JV. (93) 
Here 

A N 

S(r) = -J2Mr)^(w(m,r)) . (94) 
^ i=i 

Once all the initial amplitudes Ai(0) of the holes are 
known, the effective dynamical system (|93[) provides a 
complete description of the problem. The conservation 
law (j2"2"]l of the original Eq. (fT^|) becomes an integral of 
motion of the dynamical system (I93[) : 

N 

$> 2 (t) = t = const. (95) 

i—l 

Equations (|9"3"|) - (f9!>]) are similar to (the discrete version 
of) the Lifshitz-Slyozov theory of Ostwald ripening [35j ] . 
and their properties give a qualitative explanation to the 
properties of coarsening observed in Fig. 1101 Indeed, the 
holes with amplitudes greater than the (time-dependent) 
critical amplitude A cr {r) — (3/2) S'~ 1 (r) grow in the am- 
plitude, while holes with amplitudes less than A cr (r) de- 
crease their amplitude and disappear. As A cr (r) grows 
with time, the holes that previously grew in the ampli- 
tude begin to decrease their amplitude and finally disap- 
pear. 
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A natural further step is to assume N 1, treat the 
hole amplitude as a continuous variable and deal with the 
probability distribution F(A, r) of the hole amplitudes A 
at time r. The corresponding theory can be formulated 
in the spirit of the Lifshitz-Slyozov theory of Ostwald 
ripening, and we present it in Appendix B. How does 
this theory compare with numerical simulations? Fig- 
ure [TT] presents some quantitative characterization of the 
hole coarsening dynamics for the numerical simulation 
shown in Fig. [TTJ] Shown are the time histories of (w) 
(panel a), of the total number of holes in the system N 
(panel b) and of the sum of the hole amplitudes squared 
(panel c) for the simulation shown in Fig. [TDJ [Because 
of the noisy initial condition, it takes some time for well- 
defined holes to nucleate. We started the hole count at 
the time when the total number of the local maxima of 
w(m) became equal, for the first time, to the total num- 
ber of m-intervals where w was less than a prescribed 
small threshold 1CP 4 .] One can immediately see on the 
lower panel of Fig. [11] that the conservation law fM]) is 
not obeyed in this simulation. It is not surprising, there- 
fore, that other quantitative predictions of our Lifshitz- 
Slyozov-type theory, see Appendix B, are also not sup- 
ported by this simulation. Most directly, the shape of an 
individual hole does not agree with that assumed in the 
ansatz The holes observed in this "generic" sim- 

ulation have a more complicated structure, and are not 
characterizable by a single parameter such as Aj(r). 




120r 

80 
40 



250 500 7501000 



0. 



250 500 7501000 
x 



3x10 
< 2x10' 
1xl0 5 



250 500 7501000 



FIG. 11: The time histories of (w) (a), the number of holes 
TV (b) and the sum of the hole amplitudes squared (c) for the 
"generic" simulation (starting from a small amplitude noise) 
shown in Fig. [TUJ 
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FIG. 12: Coarsening of holes when starting from the ansatz 
(f92|) with No = 2 x 10 4 holes. The initial hole amplitudes 
Ai are randomly distributed according to a (positive) half- 
gaussian with variance 1. This distribution is normalized by 
the condition A^ = 3£/8. Shown is a small fragment of 
the system of rescaled length/mass C = 10 6 . 



It is therefore remarkable, that the ansatz does 
describe a stable regime of coarsening. That is, if one 
starts the simulation, at r = 0, with an ensemble of holes 
with different amplitudes, describable by the ansatz (|92p . 
the ansatz continues to hold and, moreover, the system 
approaches the simple scaling regime predicted by our 
theory of Lifshitz-Slyozov type. The results of one such 
simulation are presented in Figs. [T3]- [TJ] Here the holes 
were placed at a (sufficiently large) equal distance from 
each other, and the initial hole amplitudes Ai were cho- 
sen randomly from a positive half-gaussian with variance 
1. One can see a hole coarsening process in Fig.[T2J holes 
with a larger amplitude (that is, with less gas) grow (that 
is, loose gas) at the expense of holes with a smaller am- 
plitude. The time histories of (w) and the number of 
holes N, presented in Fig.[T31 resemble those for the pre- 
viously described "generic" simulation. The behavior of 
the sum Af is, however, dramatically different: here 
the conservation law ([§5]) is obeyed with a 1 percent accu- 
racy. A closer inspection of the time histories of (w) and 
N(t) (see Fig. [FT| shows that, at late times, these quanti- 
ties agree with the theoretical predictions from Eqs. (|B5|) 
(with /«i = 1) and (|B10|) . presented in Appendix B. In- 
deed, by using only one adjustable parameter: the time 
shift Tf, related to the time of approaching the scaling 
regime, we obtained good agreement for the two different 
quantities. We also checked (not shown) that, at differ- 
ent times, the shapes of individual holes are very well 



described by the cosh 

Jai. 



profile assumed in the ansatz 



VIII. SUMMARY AND DISCUSSION 

We have developed a nonlinear theory of low Mach 
number channel flows of freely cooling dilute granular 
gases with nearly elastic particle collisions. We focused 
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FIG. 14: A comparison of the time histories of (w) and N 
from Fig. [13] with theoretical predictions. The solid line 
in panel a shows the numerical results for l/{w) versus time 
r. The dashed line shows our theoretical prediction for late 
times: (w) = l/(r — Tf), where Tf is an adjustable param- 
eter (in this simulation Tf ~ 23.9). Plotted in panel b is 
the numerical result for iV -1 / 2 versus r (the solid line), and 
the theoretical prediction AT _1 ' 2 (r) = s/3/C (r — r/) with 
no additional adjustable parameters (the dashed line). Here 
V3/Z ~ 1.73 x 10~ 3 . The noise, evident in panel b at late 
times, is due to a small number of holes at those times. 



on the case when the sound travel time through the sys- 
tem is much shorter than the cooling time and the heat 
diffusion time. Then, after a brief transient, the gas pres- 
sure becomes (almost) uniform in space. This makes 
it possible to reduce the granular hydrodynamic equa- 
tions, in Lagrangian coordinates, to a single nonlinear 
and nonlocal equation of a reaction-diffusion type. With 
heat diffusion neglected, the reduced equation becomes 
integrable, and any inhomogeneous initial condition pro- 
duces a finite-time density blowup. The density blowup 
has the same universal features at singularity as those 
exhibited by a family of exact solutions of the full set 
of ideal hydrodynamic equations [ij], HI]- The density 
blowup, however, is arrested by the heat diffusion. As 
a result a novel, inhomogeneous cooling state (ICS) of 
the gas emerges which has a time-independent density 
profile. For channels of an intermediate length that we 
considered, the ICS represents a global attractor of the 
system. Both its structure, and the late-time relaxation 
towards it are determined by a single dimensionless pa- 
rameter C which is of the order of the ratio of the channel 
length to the critical length predicted by the linear the- 
ory of instability of the homogeneous cooling state. The 
energy decay of the ICS differs considerably from Haff 's 



law: the characteristic decay time diverges with the size 
of the system as C 1 ^ 2 , see Eq. (|78|) . At large £, the max- 
imum density of the ICS grows exponentially with C 
Therefore, for sufficiently long channels (the rest of pa- 
rameters being fixed), the dilute gas assumption breaks 
down, and close packed regions emerge. 

For C » 1 the cooling dynamics proceeds as a compe- 
tition between "holes" . This competition is quite similar 
to Ostwald ripening. In the simple case when the initial 
state consists of N well separated holes ~ cosh~ 2 (m/2), 
the analogy with Ostwald ripening becomes complete, 
as the "hole ripening" statistics exhibits a simple dy- 
namic scaling behavior and is describable by a variant of 
the Lifshitz-Slyozov theory. Here, in analogy with other 
phase ordering systems with a conserved order parame- 
ter, the probability distribution of the holes with respect 
to their amplitudes approaches, at long times, the spe- 
cial (limiting) self-similar solution, that is analytic at the 
edge of its (compact) support. However, for a generic, 
noisy initial condition, the competing holes have a more 
complicated structure than that described by the ansatz 
(f92|) . This brings about a lack of simple dynamic scaling. 
A theory of this regime has yet to be developed. 

In the light of the above results, a non-linear devel- 
opment of the clustering instability of the HCS, for in- 
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termediate channel lengths, is but a particular case of 
our low Mach number theory. Ultimately, the instabil- 
ity transforms an (almost) homogeneous initial gas den- 
sity profile into an inhomogeneous but time-independent 
density profile: the ICS described above. For L ^> 1 
this transformation occurs through an intermediate state 
with many holes (and many clusters). 

It would be interesting to investigate the ICSs, and 
relaxation toward them, in MD simulations. To directly 
test our low Mach number theory, one should choose the 
MD simulation parameters so as to guarantee the length 
scale separation l cr < L <C l s assumed here. We stress 
that this hierarchy of length scales demands nearly elastic 
particle collisions: \f\ — r 2 <C 1. In addition, the channel 
length L should not be too large so that the theoretically 
predicted maximum gas density in the ICs is still small 
compared to the close packing density of spheres. 

It is worth noticing that, in all asymptotic cooling 
regimes of an inhomogeneous gas that we have investi- 
gated, the energy decays slower than in the case of a 
HCS. Haff's cooling law, therefore, provides an upper 
bound on the energy decay rate. In fact, this is a general 
theorem, universally valid for a low Mach number flow. 
Indeed, according to Eq. (fT2")l . the logarithmic derivative 
of the pressure (and, therefore, of the total energy) is pro- 
portional to —(w). For a HCS (w) — 1, whereas for any 
ICS (w) < 1, by virtue of the Cauchy-Schwarz inequality 
and the identity (w 2 ) = 1. 

What can be said about the opposite, long-wavelength 
limit, A > ! s , where A is the characteristic length scale 
of the initial perturbations? Alth oug h there has been 
some progress in this case [lU, [TH, [17] . a complete un- 
derstanding of the dynamics and structure of the flow is 
still lacking. It should be possible to derive a different 
reduced model in that limit, and see whether the popular 
"pressure instability scenario" [H is at work there. (It is 
clear that the pressure instability scenario is irrelevant 
in the intermediate wavelength limit, considered in the 
present paper.) 

Note that the ICSs, that we have discovered here, are 
exact solutions of the full set of granular hydrodynamic 
equations ©-((4]) for a nearly elastic dilute gas, without 
any reductions. Therefore, a question arises on whether 
the ICS represents an attractor in the general case, in- 
cluding the long wavelength limit. A complete (unre- 
duced) linear stability analysis around the "hole" asymp- 
totics (jnOJ) could be the first step in an attempt to answer 
this question. Such an analysis can be complemented by 
numerical hydrodynamic simulations of nonlinear cooling 
flows, so as to elucidate possible effects of shock waves 
on the (nonlinear) stability of the ICS. 

Does this work, limited to channel flows, have any rel- 
evance to the shearing/clustering instability of a freely 
cooling granular gas in fully multi-dimensional geome- 
tries? To begin with, the low Mach number theory can 
be extended to the higher dimensions, once the char- 
acteristic sound travel distance l s is much larger than 
all system dimensions. This extension should take into 



a proper account the flow vorticity, in much the same 
way as it was done in Ref. [2?J where a two-dimensional 
low Mach number flow of an ideal gas, driven by the 
heat diffusion, was investigated. Although not very sim- 
ple, such a reduced description (with the acoustic modes 
eliminated) will be advantageous compared to the full set 
of multi-dimensional hydrodynamic equations. Further- 
more, the novel ICSs of the granular gas (that represent 
exact solutions of the unreduced granular hydrodynamic 
equations) may have multi-dimensional analogs. Find- 
ing these analogs, and investigating their stability with 
respect to multi-dimensional perturbations which have 
both potential, and solenoidal velocity components, can 
be a natural next step in developing a more complete 
nonlinear theory of the shearing/clustering instability. 
The channel flow theory developed here (see also Refs. 
[H, [H, H3, [H, H3| ) sets the ground for the future work. 
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Appendix A. Numerical scheme 

We employed the following implicit finite difference 
scheme for a numerical solution of Eq. (11!: 
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(Al) 



where St is the time step, Wi = w{rrii,T + <5r), u)i = 
w{rrii,T). A standard discretization Dwi of the diffusion 
term was used: for the periodic BCs we put 
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h 2 ' 
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w i — 2 w n + w n _ i 



i=l, 

1 < i < n . 
i = n . 



where h = C/n is the grid size. The approximation error 
of this scheme is 0(8t 2 ) in St and 0(h 3 ) in h. Note that 
the scheme conserves exactly the discrete version of the 
conservation law ([22jl . (wi(T) 2 ) = (1/n) Y17=i w i = 1> 
once (wi(0) 2 ) = 1. 

We solved the set of nonlinear algebraic equations (jAlj) 
by an iteration procedure based on Newton's method. To 
obtain, after linearization, a standard cyclic tridiagonal 
system, we used the values of Wi, entering the sum 52 . Wi, 
from the previous iteration. We demanded that the resid- 
ual (the maximum of the absolute value of the difference 
between the left and right hand sides of the equations 
after the iteration process) be less then 10 -13 . Because 
of the finite residual, this procedure conserved the mean 
square of w with an almost machine precision, but not 
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exactly. Therefore, we enforced an even stricter conser- 
vation by adding, at each time step, a constant c to the 
numerical solution Wi found with the iteration procedure. 
The value of c is determined as follows. We represent the 
(yet unknown) corrected values Wi as u>i — Wi + c. Then 



1 



w 2 > = -£ 

11 ' 



(w\ + 2cwi + (?) = (w 2 ) + 2c(w) + c 2 



Now we demand that the right hand side be equal to 1. 
Neglecting the c 2 term, we find 



1 



(w 



2{w) 

We always obtained |c| < 10 -14 in our computations. 
This justifies neglecting the c 2 term. 

The typical set of parameters for the investigation of 
relaxation towards a stationary single hole asymptotics 
(Hi]) was C = 50 and n = 2.5 x 10 4 , so h = 2 x 1(T 3 . 
In the hole coarsening simulations we used £ = 10 6 and 
n = 2.8 x 10 6 , so h ~ 0.36. In all cases the time step was 
chosen to be St = h 2 . 



Appendix B. Hole coarsening in the spirit of the 
Lifshitz-Slyozov theory 

Here we treat the hole amplitude (see Section IVTl]) as a 
continuous variable and deal with the probability distri- 
bution F(A, r) of the hole amplitudes A at time r. The 
total number of holes N(t) = F(A,r)dA > 1. As 
there is no nucleation of new holes and no hole mergers, 
the evolution of F(A, r) is described, in the spirit of the 
Lifshitz-Slyozov theory [351 ] . by a continuity equation in 
the space of hole amplitudes: 



dF 



d_ 
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SA--)F 



= 0. 



where 



and 



S(t) = ^ AF(A,r)dA, (B2) 



/ A 2 F(A, r)dA = -^ = const . (B3) 

Equations similar to Eqs. (|B1|) - (|B3[) have appeared in the 
context of the Lifshitz-Slyozov model of Ostwald ripen- 
ing |35fl and its analogs for different transport mecha- 
nisms [25l I3&1 l37l [H, H§| . In those systems one is usually 
interested in the question of whether or not the prob- 
ability distribution F(A, r) approaches, at late times, a 
self-similar shape. A simple power counting in Eqs. (|Bip ~ 
(|B3t yields 



(B4) 



where <&{rj) > is the (yet unknown) rescaled distribu- 
tion, and the coefficient C/4 is chosen for convenience. 
Using Eqs. JB2J) and JB3J, we obtain 
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(B5) 



respectively. Here fik is the fc-th moment of the rescaled 
distribution: /ij. = J °° rj k ^{rj)drj. One can already see 
that the total number of holes N{t) goes down as t~ 2 , 
while both the average hole amplitude A(t) and the crit- 
ical amplitude A cr (r) grow linearly with r. The pre- 
factors of these power laws will be determined once $(77) 
is found. Plugging Eq. (|B4[1 and the first of Eqs. (|B5D 
into Eq. (|B1[) we obtain an ordinary differential equation 
for $(77): 



(Mi 



3)$ = 0. 



(B6) 



whose solution is elementary. As in other variants of 
the LS-theory, we obtain here a whole family of shape 
functions $^(77), parameterized by the first moment 
The solutions exist, with finite moments, for 1 < fii < 00. 
For ui > 1 the solutions have finite support: 
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if V > Vm , 



(B7) 

where r\ m — (3/2) (/ii — 1) 1 . The constant B fil can be 
determined from the second of Eqs. (|B5|) (that plays the 
role of a normalization condition): 



B„ 



2^1 _ Ml + 1 

2^1-1 3 Mi-i fii(fii 



(Bl) N(t) 



This yields, for [i\ > 1, 
£/xi(l + fix) 



6r 2 



A(r) 



(^i-i) 2 r(^) 



A ct {t) = (3r)/(2 Ml ). 



(B8) 



For 1 < Hi < 3, the solutions (|B7[) vanish at 77 = i] m , 
whereas for ui > 3 they diverge at 7/ = r\ m . As all the 
moments /j^- remain finite, the diverging distributions are 
legitimate. 

For /ii = 1 we obtain a limiting solution $1(77) = 
(16/9) exp(—4?7/3) that has an infinite support < 77 < 
00. The self-similar probability distribution (|B4p be- 
comes 



4£ 

F ( A ' r )= 9^3 ex P 



AA 
' 3r 



(B9) 



In this case 



£ T/ , 3r 



3r 



N{t) = — , A(t) = T , and A cr (r) = — . (BIO) 



These expressions also follow from Eqs. (|B8[) in the limit 
of hi —> 1. 
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FIG. 15: Rescaled distributions of the hole amplitudes <I> M1 (77) 
for fi\ — 2.5, 2, 1.5, and 1 (the latter corresponding to the 
limiting distribution). Smaller /ii are shown by thicker lines. 



Figure [15] depict the rescaled distributions <j> Ml (77) for 
four values of the parameter fix. Selection of the "cor- 
rect" self-similar solution out of the family of solution 
represents a subtle problem that was resolved only re- 
cently. It turns out that the selection is only made by 
(a certain feature of) the initial condition F(A, r = 0) 
[H, M, H [H HI. If F(A,t = 0) has compact sup- 
port, the similarity solution, if any, is selected by the 
behavior of F(A, t — 0) near the supremum A max of the 
support. If F(A, t — 0) has a power-law asymptote near 
A m ax, the exponent of this power law selects one of the 



solutions from the family (|B7|) . If F(A,t — 0) goes to 
zero exponentially fast at A — > A max (or if the support 
of F(A,t = 0) is infinite), the limiting solution (|B9[) is 
selected. 

This sensitivity to initial conditions shows a certain 
lack of robustness of the Lifshitz-Slyozov model and its 
analogs like our Eqs. (|B1[) - (|B3[) . As a remedy, one has 
to account for an additional physics (that may be less 
universal and more system-dependent). For example, 
in the context of the interface-controlled Ostwald ripen- 
ing strong selection is achieved via an account of direct 
droplet merger events [ill ]. 

As we show in Section IVII1 the Lifshitz-Slyozov-type 
model does not agree with numerical simulations that 



start from generic initial conditions. However, if one 
starts the simulation with an assembly of holes, describ- 
able by the ansatz (j9"2")l , the ansatz continues to hold, 
and the system approaches the simple scaling regime pre- 
dicted by the Lifshitz-Slyozov-type theory. Therefore, we 
want to pursue the ansatz (|92p a bit further, as it pro- 
vides an interesting, though non-generic, characterization 
of the hole coarsening. We assume that the limiting dis- 
tribution (|B9[) . corresponding to fix = 1, is selected and 
use Eq. (j2"0|) and the relation S(t) = 1/t to find the cor- 
responding scaling behavior of the gas pressure p(r). We 
obtain 



l dp 
p(r) dr 



.2 7 S( r ) = -_L 



2 _2 

r 



(Bll) 



which yields p(r) — po(ro/r) 27 , where To is an effective 
"initial" time, and po = p{tq). Using Eq. (|2l"j) . we find the 
following relation between the original (physical) time t 
and the new time t: 



27 r 



7+1 



(7 



As a result, 

m = 

Po 



27 TO 



( 7 +l)A/#%i 



■2-1 

■y+l 



2-1 
1 + 1 



where to = t(ro). Now, in the low Mach number regime 
we have been dealing with throughout this paper, the 
total energy of the gas decays in (almost) the same 

2~/ 

way as the pressure, so E to t(t) ~ t 1+1 ■ We obtain 
Etotit) ~ r 4 / 3 and E tot (t) ~ i~ 5 / 4 in 2d (disks) and 
3d (spheres), respectively. Again, the cooling dynamics 
proceeds slower than that predicted by Haff's law (|25[) . 
We checked that the same conclusion holds for any fix, 
that is for all possible self-similar distributions of the hole 
amplitudes. 
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